Maximizing as minimum
Adapted from: (Floudas et al., 1999; Section 4.10), (Laurent, 2008; Example 6.23) and (Lasserre, 2009; Table 5.1)
Introduction
Consider the polynomial optimization problem from (Floudas et al., 1999; Section 4.10) of minimizing the linear function $-x_1 - x_2$ over the basic semialgebraic set defined by the inequalities $x_2 \le 2x_1^4 - 8x_1^3 + 8x_1^2 + 2$, $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$,
using DynamicPolynomials
@polyvar x[1:2]
p = -sum(x)
using SumOfSquares
f1 = 2x[1]^4 - 8x[1]^3 + 8x[1]^2 + 2
f2 = 4x[1]^4 - 32x[1]^3 + 88x[1]^2 - 96x[1] + 36
K = @set x[1] >= 0 && x[1] <= 3 && x[2] >= 0 && x[2] <= 4 && x[2] <= f1 && x[2] <= f2Basic semialgebraic Set defined by no equality
6 inequalities
x[1] ≥ 0
3 - x[1] ≥ 0
x[2] ≥ 0
4 - x[2] ≥ 0
2 - x[2] + 8*x[1]^2 - 8*x[1]^3 + 2*x[1]^4 ≥ 0
36 - x[2] - 96*x[1] + 88*x[1]^2 - 32*x[1]^3 + 4*x[1]^4 ≥ 0
As we can observe below, the bounds on x[2] could be dropped and optimization problem is equivalent to the maximization of min(f1, f2) between 0 and 3.
xs = range(0, stop = 3, length = 100)
using Plots
plot(xs, f1.(xs), label = "f1")
plot!(xs, f2.(xs), label = "f2")
plot!(xs, 4 * ones(length(xs)), label = nothing)
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy gives a lower bound of -7`
model4 = solve(4)-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 31
constraints = 40
nnz(P) = 0
nnz(A) = 73
cones (total) = 6
: Zero = 1, numel = 10
: PSDTriangle = 5, numel = (6,6,6,6,6)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 1.2590e-01 1.2590e-01 8.33e-17 7.72e-01 5.10e-01 1.00e+00 1.51e+00 ------
1 4.8332e-01 5.0333e-01 2.00e-02 1.61e-01 9.82e-02 1.86e-01 3.32e-01 8.02e-01
2 1.3995e+00 1.5138e+00 8.17e-02 8.80e-02 3.22e-02 2.10e-01 1.08e-01 7.85e-01
3 2.7854e+00 2.9400e+00 5.55e-02 1.81e-02 4.48e-03 1.83e-01 1.58e-02 8.78e-01
4 4.0242e+00 4.2592e+00 5.84e-02 1.24e-02 1.85e-03 2.61e-01 7.20e-03 7.46e-01
5 5.4090e+00 5.5096e+00 1.86e-02 3.46e-03 4.31e-04 1.09e-01 1.81e-03 8.27e-01
6 6.2811e+00 6.3385e+00 9.13e-03 1.45e-03 1.50e-04 6.15e-02 6.79e-04 7.22e-01
7 6.9685e+00 6.9720e+00 5.12e-04 5.52e-05 5.40e-06 3.75e-03 2.52e-05 9.69e-01
8 6.9993e+00 6.9994e+00 1.29e-05 1.39e-06 1.35e-07 9.46e-05 6.33e-07 9.75e-01
9 7.0000e+00 7.0000e+00 3.48e-07 3.79e-08 3.68e-09 2.56e-06 1.72e-08 9.73e-01
10 7.0000e+00 7.0000e+00 1.18e-08 1.30e-09 1.25e-10 8.65e-08 5.85e-10 9.66e-01
11 7.0000e+00 7.0000e+00 2.32e-10 1.90e-10 3.31e-12 1.73e-09 1.18e-11 9.74e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 2.34ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -7.00000e+00
Dual objective value : -7.00000e+00
* Work counters
Solve time (sec) : 2.33684e-03
Barrier iterations : 11The second level improves the lower bound
model5 = solve(5)-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 108
constraints = 128
nnz(P) = 0
nnz(A) = 266
cones (total) = 7
: Zero = 1, numel = 21
: Nonnegative = 1, numel = 2
: PSDTriangle = 5, numel = (21,21,21,21,21)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 2.8313e-01 2.8313e-01 0.00e+00 8.12e-01 7.55e-01 1.00e+00 1.82e+00 ------
1 5.4905e-01 5.6026e-01 1.12e-02 1.47e-01 1.42e-01 1.86e-01 4.16e-01 8.05e-01
2 1.0075e+00 1.0482e+00 4.04e-02 3.92e-02 4.05e-02 8.87e-02 1.08e-01 8.48e-01
3 1.7352e+00 1.7929e+00 3.33e-02 8.93e-03 6.29e-03 7.10e-02 1.55e-02 9.20e-01
4 2.2703e+00 2.3425e+00 3.18e-02 3.95e-03 1.67e-03 7.97e-02 4.20e-03 8.20e-01
5 2.5204e+00 2.5826e+00 2.47e-02 2.72e-03 9.42e-04 6.81e-02 2.42e-03 6.21e-01
6 3.3688e+00 3.4392e+00 2.09e-02 7.06e-04 1.40e-04 7.25e-02 3.73e-04 8.78e-01
7 4.1613e+00 4.2336e+00 1.74e-02 3.97e-04 3.95e-05 7.36e-02 1.14e-04 8.80e-01
8 4.9013e+00 4.9482e+00 9.56e-03 1.79e-04 1.33e-05 4.76e-02 4.21e-05 8.78e-01
9 5.7265e+00 5.7548e+00 4.95e-03 6.71e-05 3.22e-06 2.86e-02 1.12e-05 8.50e-01
10 5.7843e+00 5.8135e+00 5.04e-03 6.48e-05 2.90e-06 2.95e-02 1.03e-05 2.28e-01
11 6.4423e+00 6.4536e+00 1.75e-03 1.41e-05 5.51e-07 1.13e-02 2.07e-06 8.34e-01
12 6.6492e+00 6.6500e+00 1.19e-04 8.39e-07 3.36e-08 7.94e-04 1.26e-07 9.69e-01
13 6.6662e+00 6.6663e+00 3.10e-06 2.18e-08 8.74e-10 2.08e-05 3.28e-09 9.74e-01
14 6.6667e+00 6.6667e+00 8.84e-08 6.27e-10 2.49e-11 5.93e-07 9.36e-11 9.72e-01
15 6.6667e+00 6.6667e+00 2.50e-09 5.37e-10 1.67e-12 1.69e-08 2.84e-12 9.37e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 7.01ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -6.66667e+00
Dual objective value : -6.66667e+00
* Work counters
Solve time (sec) : 7.00918e-03
Barrier iterations : 15The third level finds the optimal objective value as lower bound...
model7 = solve(7)-------------------------------------------------------------
Clarabel.jl v0.10.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 288
constraints = 323
nnz(P) = 0
nnz(A) = 739
cones (total) = 8
: Zero = 1, numel = 36
: PSDTriangle = 7, numel = (55,55,55,55,...,55)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 3.3234e-01 3.3234e-01 2.22e-16 8.45e-01 6.51e-01 1.00e+00 1.51e+00 ------
1 4.7370e-01 4.8104e-01 7.34e-03 1.45e-01 1.49e-01 1.77e-01 3.55e-01 8.03e-01
2 7.3313e-01 7.6805e-01 3.49e-02 4.83e-02 6.81e-02 1.12e-01 1.54e-01 7.48e-01
3 1.1439e+00 1.1684e+00 2.14e-02 8.88e-03 1.34e-02 4.09e-02 2.71e-02 8.65e-01
4 1.5684e+00 1.6151e+00 2.98e-02 3.54e-03 4.11e-03 5.81e-02 7.91e-03 8.51e-01
5 1.8613e+00 1.8904e+00 1.57e-02 1.63e-03 1.58e-03 3.50e-02 3.10e-03 6.74e-01
6 2.1699e+00 2.1913e+00 9.85e-03 5.93e-04 4.34e-04 2.40e-02 8.88e-04 8.03e-01
7 2.5021e+00 2.5220e+00 7.93e-03 2.14e-04 1.12e-04 2.11e-02 2.37e-04 8.01e-01
8 2.7896e+00 2.8141e+00 8.76e-03 1.35e-04 4.21e-05 2.54e-02 9.20e-05 8.08e-01
9 3.2619e+00 3.2812e+00 5.91e-03 4.58e-05 9.63e-06 1.96e-02 2.16e-05 7.88e-01
10 3.6591e+00 3.6805e+00 5.86e-03 3.22e-05 3.07e-06 2.17e-02 7.16e-06 9.60e-01
11 4.0143e+00 4.0278e+00 3.36e-03 1.51e-05 1.27e-06 1.36e-02 3.07e-06 7.92e-01
12 4.3776e+00 4.3897e+00 2.76e-03 8.80e-06 5.15e-07 1.22e-02 1.31e-06 6.54e-01
13 4.8832e+00 4.8912e+00 1.64e-03 2.91e-06 1.15e-07 8.04e-03 3.30e-07 8.39e-01
14 5.0180e+00 5.0266e+00 1.70e-03 2.39e-06 7.30e-08 8.58e-03 2.09e-07 6.21e-01
15 5.3306e+00 5.3332e+00 4.93e-04 6.71e-07 2.20e-08 2.64e-03 6.20e-08 9.90e-01
16 5.4834e+00 5.4837e+00 5.85e-05 8.85e-08 2.99e-09 3.22e-04 8.40e-09 8.73e-01
17 5.5049e+00 5.5049e+00 1.02e-05 1.10e-08 3.75e-10 5.63e-05 1.05e-09 9.61e-01
18 5.5079e+00 5.5079e+00 5.58e-07 5.98e-10 2.05e-11 3.08e-06 5.75e-11 9.46e-01
19 5.5080e+00 5.5080e+00 3.50e-08 3.43e-11 1.17e-12 1.93e-07 3.27e-12 9.88e-01
20 5.5080e+00 5.5080e+00 8.03e-10 1.42e-12 3.81e-14 4.43e-09 7.49e-14 9.77e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 30.0ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -5.50801e+00
Dual objective value : -5.50801e+00
* Work counters
Solve time (sec) : 2.99861e-02
Barrier iterations : 20...and proves it by exhibiting the minimizer.
ν7 = moment_matrix(model7[:c])
η = atomic_measure(ν7, 1e-3) # Returns nothing as the dual is not atomicAtomic measure on the variables x[1], x[2] with 1 atoms:
at [2.3295201628688242, 3.17849313493404] with weight 0.9999999696085651We can indeed verify that the objective value at x_opt is equal to the lower bound.
x_opt = η.atoms[1].center
p(x_opt)-5.508013297802864We can see visualize the solution as follows:
scatter!([x_opt[1]], [x_opt[2]], markershape = :star, label = nothing)
This page was generated using Literate.jl.